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Abstract. We investigate an XY spin-glass model in which both spins and 
interactions (or couplings) evolve in time, but with widely separated time-scales. For 
large times this model can be solved using replica theory, requiring two levels of 
replicas, one level for the spins and one for the couplings. We define the relevant order 
parameters, and derive a phase diagram in the replica-symmetric approximation, which 
exhibits two distinct spin-glass phases. The first phase is characterized by freezing of 
the spins only, whereas in the second phase both spins and couplings are frozen. A 
detailed stability analysis leads also to two distinct corresponding de Almeida-Thouless 
lines, each marking continuous replica-symmetry breaking. Numerical simulations 
support our theoretical study. 



PACS numbers: 75.10.Nr, 05.20.-y, 64.60.Cn 
1. Introduction 

The study of coupled dynamics of fast Ising spins and slow couplings has received 
considerable interest recently (see e.g. |J-0 and references therein), stimulated by 
considerations of simultaneous learning and retrieval in recurrent neural networks and 
the influence of slow atomic diffusion processes in disordered magnetic systems. 

Generalizing spin systems by taking their interactions to be (slowly) time dependent 
was first considered in 0, as a mechanism with which to restore broken ergodicity at 
low temperature in the SK model |7|]. Another conceptually similar process, but now 
describing slow and deterministic synaptic modification in neural systems, driven by 
averages over neuron states, was first introduced in ||. Explicit stochastic dynamical 
laws for the interactions were defined in ||, where spin-glass models with coupled 
dynamics were studied within replica mean-field theory. It turned out that the replica 
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dimension in such models has a direct physical interpretation as the ratio of two 
temperatures characterizing the stochasticity in the spin dynamics and the coupling 
dynamics, respectively. Later it was shown that the case of negative replica dimension 
represents an over-frustrated system 0. In a similar spirit, neural network models 
with a coupled dynamics of fast neurons and slow neuronal connections were treated in 
0-0- 

In this paper the results previously obtained by others for Ising spin models are 
further extended to a classical XY spin glass with dynamic couplings, whose continuous 
spin variables are physically more realistic than Ising ones. Moreover, the XY model is 
closely related to models of coupled oscillators |L5], of which the neural network version 
fl6| provides a phenomenological description of neuronal firing synchronization in brain 
tissue. In particular, we examine the effects of including an explicit frozen randomness 
into the dynamics of the interaction weights. 

The model is solved using the replica formalism. Relevant order parameters are 
defined and a phase diagram is obtained upon making the replica-symmetric Ansatz. 
Similarly to the Ising case, we find two different spin glass phases in addition to a 
paramagnetic phase. One spin-glass phase exhibits freezing of the spins in random 
directions, but on the time-scale of the coupling dynamics these 'frozen directions' still 
continue to change. A second spin-glass phase exhibits freezing of the spins as well 
as of the couplings, such that even on the large time-scales the 'frozen directions' of 
the spins remain stationary. We perform a detailed stability analysis and calculate 
the de Almeida-Thouless (AT) lines |DJ (of which here there are two types), where 
continuous transitions occur to phases of broken replica symmetry. A brief preliminary 
account of the first part of the present work has been presented in ||18|| . Finally, we 
discuss and tackle the problem of simulating this model numerically. 

The remainder of this paper is organized as follows. In section 2 the classical XY 
spin glass model with coupled dynamics is defined. In section 3 the order parameters 
are calculated in the replica symmetric (RS) Ansatz and a phase diagram is presented. 
As is well known, the solutions of this Ansatz are not always stable against replica 
symmetry breaking (RSB). Therefore, the lines of instability are calculated in section 4. 
Finally, section 5 presents the results of the numerical simultions of this model, followed 
by a concluding discussion in section 6. The appendix describes all eigenvectors and 
eigenvalues of the Hessian matrix determining the stability of the replica symmetric 
solutions. 



2. The model 



We consider a system of N classical two-component spin variables Si = (cos 0$, sin 0$), 
i = 1...N, with symmetric couplings (or exchange interactions) J^-, taken to be of 
infinite range. In contrast to the standard XY spin glass, these couplings are not static 
but are allowed to evolve in time, albeit slowly. The spins are taken to have a stochastic 
Glauber-type dynamics such that for stationary choices of the couplings the microscopic 
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spin probability density would evolve towards a Boltzmann distribution 

P({5,}, {Jij}) ~ exp[-/3tf ({£,}, {Jij})} (1) 

with the standard Hamiltonian 

H({S i },{J ij }) = - y £JkiS k -S t (2) 

and with inverse temperature (3 = T -1 , where k, I G {1, . . . , N}, and where, at least for 
the purpose of the dynamics of the spins, the {Jij} are to be considered as quenched 
variables. 

We remark that this system is equivalent to a system of N coupled oscillators with 
phases 0j |T3|, whose time evolution is described by a Langevin equation 

j t e t = Ja MOj - Oi) + JjUt) , (3) 

where the £i(t) are defined as independent white noise variables, drawn from a Gaussian 
probability distribution with 

(&(*)> = o (&(t)&(0> = M(*-f)- (4) 

In our model, the couplings also evolve in a stochastic manner, partially in response 
to the states of the spins and to externally imposed biases. However, we assume that 
the spin dynamics is very fast compared to that of the couplings, such that on the time- 
scales of the couplings the spins are effectively in equilibrium (i.e. we take the adiabatic 
limit). For the dynamics of the couplings the following Langevin form is proposed : 

p^ {S,-S J ) + K,, _^ + ^l . ,... A . (5) 

The term (Si ■ Sj), representing local spin correlations associated with the coupling J^, 
is a thermodynamic average over the Boltzmann distribution (|T]) of the spins, given 
the instantaneous couplings {Jke}- No other spins are involved, in order to retain the 
local character of the couplings. We remark that only the thermal averages (or long time 
averages) of the spin correlations play a role, rather than the instantaneous correlations, 
since the dynamics of the couplings is (by definition) sufficiently slow. External biases 
Kij = fiNBij serve to steer the weights to some preferred values. The Bij are chosen 
to be quenched random variables, drawn independently from a Gaussian probability 
distribution with mean Bq/N and variance B/N: 

(Bu-Bq/N)"' 



P( B ij) = i ~ ex P 
'2-kB/N 



2B/N 



(6) 



and are thus reminiscent of the couplings in the original SK model [[?]. Here, in 
contrast, the {B^} generate frozen disorder in the dynamics of the couplings. The 
decay term fiJij in (^) is added in order to limit the magnitude of the couplings. 
Finally, the terms r)ij(t) represent Gaussian white noise contributions, of zero mean 
and covariance (Vij(t)Vki(t')) = 2T 8i^8ji8(t — t'), with associated temperature T = 
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Appropriate factors of iV are introduced in order to ensure non-trivial behaviour in the 
thermodynamic limit iV — ► oo. 

The model exhibits three independent global symmetries, which can be expressed 
efficiently in terms of the Pauli spin matrices a x and a z : 

inversion of both spin axes : Si —Si for all i 

inversion of one spin axis : Si — > o z S{ for all % (7) 

permutation of spin axes : Si — > o~ x Si for all i . 

Upon using algebraic relations such as a x a z a x = —o z and o z a x a z = —o x we see that in 
the high T (ergodic) regime these three global symmetries generate the following local 
identities, respectively: 

(Si) = 0, (S i -a x S j ) = 0, (S i -a z S J ) = 0. (8) 

We note that the stochastic equation (|) for the couplings is conservative, i.e. it 
can be written as 

with the following effective Hamiltionian for the couplings: 

H({Jij}) = -\ log Z p {{J l3 }) + \n J2 Jit - M# £ B u Ju ■ (10) 

The term Zp{{Jif\) = Trrg, exp\flJ2k<e Jkt.Sk 4 Se] i n this expression is the partition 
function of the XY spins with instantaneous couplings {Jij}- It follows from @ that 
the stationary probability density for the couplings is also of a Boltzmann form, with 
the Hamiltonian (|1~0|), and that the thermodynamics of the slow system (the couplings) 
are generated by the partition function = fTik<edJki exp[— $H({ Jij})], leading to 
(modulo irrelevant prefactors): 



%p = J n i z p({ j ij})] n ex p 



k<i 1 k<i 



k<e 

In contrast to the more conventional spin systems with frozen disorder, where the 
replica dimension n is a dummy variable, here we find that n is given by the ratio 
n = (3/(3, and can take any real non-negative value. The limit n — > corresponds to a 
situation in which the coupling dynamics is driven purely by the Gaussian white noise, 
rather than by the spin correlations. Therefore, in this limit the model is equivalent 



to the XY model with stationary couplings formulated, as in [I9|. For n = 1 the two 
characteristic temperatures are the same, and the theory reduces to that corresponding 
to the exchange interactions being annealed variables. In the limit n — » oo the influence 
of spin correlations on the coupling dynamics dominates, and the couplings J^ only 
fluctuate modestly (if at all) around their mean values ((Si ■ Sj) + K i3 )/ [iN. 
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3. Statics 



We define the disorder-averaged free energy per site 

/ = -^<logi^> B , (12) 

in which (■)# denotes an average over the {-B^}- We carry out this average using the 
identity logZ^ = lim^o r -1 [iT~— 1], evaluating the latter by analytic continuation from 
integer r. Our system, characterized by the partition function Zs, is thus replicated r 
times; we label each replica by a Roman index. Each of the r functions Zz, in turn, 
is given by (|TT|) , and involves Zp({Jij}) n which is replaced by the product of n further 
replicas, labeled by Greek indices. For non-integer n, again analytic continuation is 
made from integer n. Therefore, performing the disorder average in / boils down to 
performing the disorder average of [Zg] r , involving nr coupled replicas of the original 
system: {Si} — > {Sf a }, with a = 1 . . . n and a = 1 . . . r. We obtain 

N I 1 / 2 ' 

m r )B = i n {dBn P (B iS )} m\udj% - 



Kj 



2nJ 



x Tr {Sr a } ex P 



E EK) 2 + t E E s i3 j« + p E E E b«sz ■ s% 



(13) 



where J = 1/ fi/3, and with the Gaussian probability distribution of the external biases 



Bij as given by eq. 



The Roman indices (a, b, 



run from 1 to r; the Greek ones 



(a, /3, . . .) from 1 to n. Expression ([13]) can be evaluated using the standard techniques 
of replica mean-field theory |20] . Because of the complexity of the replica structure we 
indicate the most important steps. We first perform the integrals over the couplings 
and the biases, giving 

/ \ 2 



([Z ] r ) B = Tr {5:j exp 



.Bi 



B 



/^EEE^-^ + z^EIEE^-^ 



i<j 



N 



i<j 



2N 



i<j 



(14) 



-,a(3\ 



and decouple the i- and j-components using 

ca QOt cP — ^ I ( C a P\ 

Here the quantity (S^f )j is defined as a two-dimensional unit vector with reference angle 



(15) 



equal to the difference of the reference angles of Sf a and whereas ), is defined 
as a two-dimensional unit vector with reference angle equal to the sum of both these 
angles. Upon applying the saddle-point method in the thermodynamic limit iV — > oo 
we then arrive at 

AT extr F({r<}, {a?}, {s a J }, {tf }, {tff })' 



>B 



exp 



(16) 



F({m:}, { S f}, Rf }, {tf }, {t a J}) = - l -Bf e E OS) 2 + (*: 

° a^6 a/3 



ap\2 
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afi\2 



aji/3 



a a a a 



+ logG({r<}, {sf }, {s%}, {if}, {C» 



(17) 



G({m-}, {*f }, {tf }, {C» = Tr 



,a/3 



«/3l 



+7^EE(«S-^ + *2'-r a6 



a/3 rria.fi 



aj^b a,/3 



1 , ~ 



a/3 rpa/3 



a^/3 

a a 

The parameters {m^}, {s^}, {s°f }, {if} and } introduced by this procedure 
are vectors, hence the extremum is taken over both components. They carry Greek and 
Roman replica labels. 

Those parameters which have only one Greek and Roman replica label {m®, i"), 
can be interpreted as 

1 E(M> ) R £ = &^E(^> fl - as) 



lim , 



The horizontal bar denotes thermal averaging over the coupling dynamics with fixed 
biases {Bij}. Those parameters which involve pairs of replicas either connect two distinct 
Greek replicas with a single Roman replica, (sf , if), or with two distinct Roman 
replicas, (s"f , if ). The latter vector variables can, equivalently, be expressed in terms 
of the following scalar order parameters, which measure the correlations between the 
various replicas: 



a/3 
Qab 



lim — 



3? 



u 



a/3 
ab 



N 



lim — 



Qta 



a/3 
' J ab 



lim — 



°z S ib 



B 



(20) 



At this point we remark that the order parameters wf and t>f are typical for the 
XY-model [jl9 | , and do not appear in the SK-model; comparison with (§) shows that, 
together with m° and t", they measure the breaking of the global symmetries (^). For 
simplicity we will henceforth choose B = 0. We will make the usual assumption that, 
in the absence of global symmetry-breaking forces, phase transitions can lead to at most 
local violation of the identities (§). Thus the latter will remain valid if averaged over 
all sites, at any temperature, which implies that m" = i" = and that wf = v*f, = 0. 
The spin-glass order parameters q*j. ; , on the other hand, are not related to simple global 
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symmetries, but measure the overlap of two vector spins, and serve to characterize the 
various phases. 

At this stage in the calculation we make the replica symmetry (RS) Ansatz. 
Since observables with identical Roman indices refer to system copies with identical 
couplings, whereas observables with identical Roman indices and identical Greek indices 
refer to system copies with identical couplings and identical spins, in the present 
problem the RS ansatz for the spin-glass order parameters takes the form = 
S a b {8 a /3 + qi[l — Sap}} + qo[l — S a b\. Here we remark that ■ = 1 and that, in 
the absence of global symmetry breaking forces, s"f becomes a vector of length g^f and 
reference angle 0. 

The asymptotic disorder-averaged free energy per site can now be written as 
/ = \BP 2 n 2 ql - 1 -{B + J)fn{n - l)qf -\{B + J)(3 2 n qi 



+ J Vp log | J Vq 



Tr 



5 ex P 



(3 J^Bq p-S + (3 



q-S 



(21) 



with the short-hand E = f3\] \{J +B)q\ — hBq , and where we have introduced the two- 



2\" 1 "1^ 2 

dimensional Gaussian measure 



Vp = (27r) dp x dp y exp 



(22) 



2" ' ' y ' 

The remaining two order parameters go an d q\ are determined as the solutions of the 
following coupled saddle-point equations 

2 



dx P(x) 



dx P(x) 



Jdz P{z) [IpizE)] 11 - 1 hjzE) hjzxPZ-iy/lBqo) 
Jdz P(z) [h{z~)] n I Q (zxPE-iyJ\Bti) 



Jdz P(z) [I (zE)] n - 2 [I x {z~)Yh{zxPE- l J\Bq 



Jdz P(z) [I (zZ)} n I (zxpZ^JlBq 



(23) 



(24) 



with P(x) =xe~z x 6[x], and where the functions I n (x) are the Modified Bessel functions 
of integer order [211 . 

One can give a simple physical interpretation of these order parameters in terms of 
the appropriate averages over the various dynamics 



go 



lim — y < (Si) 



B 



Qi 



lim — y ( (Si) 



B 



(25) 



It is clear that < qo < q± < 1. 

We have studied the fixed-point equations (23,24), after having first eliminated the 
parameter redundancy by putting B = 1 and J = 3. The resulting phase diagram 
in the n-T plane is shown in Fig. [I]. The appearance of two different spin-glass order 
parameters suggests that two different spin-glass phases are to be expected. Indeed, in 
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Figure 1. Phase diagram of the XY spin glass with slow dynamic couplings, drawn 
in the n-T plane with Bo = 0, B = 1 and J = 3. P: paramagnetic phase, qi = qa = 0; 
SGI: first spin-glass phase, qi > and go = (freezing on spin time-scales only); 
SG2: second spin-glass phase, qi > and qo > (freezing on all time-scales); AT lines: 
A /v. = (Roman replicon), Xq = (Greek replicon). 



addition to a paramagnetic phase (P), where go = gi = 0, we find two distinct spin-glass 
phases: SGI, where q\ > but go = 0, and SG2, where both qi > and go > 0. 
The SGI phase describes freezing of the spins on the fast time-scales only (where spin 
equilibration occurs); on the large time-scales, where coupling equilibration occurs, we 
find that, due to the slow motion of the couplings, the frozen spin directions continually 
change. In the SG2 phase, on the other hand, both spins and couplings freeze, with 
the net result that even on the large time-scales the frozen spin directions are 'pinned'. 
The SG1-SG2 transition is always second order and occurs for T = (n — l)q\ + 1/2. 
The transition SG1-P is second order for n < 2 (in which case its location is given by 
B + J = 4T 2 ), but first order for n > 2. When n further increases to n > 3.5, the 
SGI phase disappears, and the system exhibits a first order transition directly from 
P to SG2. Fig. shows for several values of n the values of the order parameters as a 
function of the temperature. 

Qualitatively, the phase diagram of the present model is very similar to that of 
the Ising spin glass with dynamic couplings 0. The main difference is the re-scaling 
by a factor two of the transition temperature from the first spin-glass phase to the 
paramagnetic phase, as has already been noticed in 

The existence of two types of spin-glass order parameters is directly related to 
the presence of quenched disorder in the couplings, which allows the latter to freeze 
in random directions at low coupling temperature T. In a model with homogeneous 
external biases [1], p^l , where no preferred direction of the couplings is assumed, one 
distinguishes (in contrast to the present situation) only the paramagnetic phase and 
the spin-glass phase SGI. Qualitatively, the transition line separating the paramagnetic 
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Figure 2. Dependence of the order parameters go (broken curve) and q\ (full curve) 
on the temperature T for various temperature ratios n, all at Bq = 0, B = 1 and 
J = 3. For n = 1.0 there is a continuous phase transition from P to SGI at T = 0.5 
and from SGI to SG2 at T = 1. For n — 2.5 the first transition occurs at T = 0.98 
while the parameter q\ drops discontinuous from 0.18 to at T = 1.03. Finally for 
n = 10 both order parameters vanish at T = 1.83 (limit value qo = 0.13, qi = 0.29) 
indicating a first order transition from P to SG2. 



phase from the spin-glass phase in the case of absent coupling disorder is the same as 
that in Fig. [l], viz. a second order transition for n < 2 given by J = 2T and a first 
order transition for n > 2. The corresponding expressions for the order parameters can 
immediately be deduced from the results above: when the quenched disorder in the 
couplings is absent, the partition function itself is self-averaging and the replica method 
is simply no longer needed. Therefore all order parameters concerning different Roman 
indices are redundant and drop out automatically, such that one ends up with only one 
spin-glass order parameter. Its explicit value is obtained by putting B = and B = 
in (0). 



4. Stability of the replica-symmetric solutions 

Additional transitions may occur in our model due to a continuous breaking of replica 
symmetry. Here we expect two distinct types of replica symmetry breaking, with respect 
to the two distinct replicas, viz. the Roman and the Greek ones. The stability of the RS 
solution is, as always, expressed in terms of the matrix of second derivatives of quadratic 
fluctuations at the saddle point [17|. We calculate all eigenvalues and their multiplicity 



following the ideas in |||, T7j- We remark that our results differ from, and improve upon 
those of H . It turns out that the (restricted) set of eigenvectors and eigenvalues given in 
satisfy only part of the relevant orthogonality conditions used in their calculations. 
In the following we present a summary of the results. More details can be found in 
Appendix A| . 



We start by rewriting (0), taking into account its invariance with respect to the 
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m<tfh tif}) = 4 bp 2 E E fe?) 2 -Ub + J)P 2 e E 



(26) 



+ log G s ({q a J},Uf}) 



G'5({g a f},{gf}) = Tr {5:} exp 



4 a/3 



(27) 



We consider small fluctuations of the order parameters around their RS saddle-point 
values 

qf = Qo + ef (a<(3) and gff = (a < 6) (28) 

and expand (^7j) up to second order in e"' 3 and r/^f. The first order terms vanish by 
construction. The coefficients of the second order terms form the so-called Hessian 
matrix and are denoted by 



d 2 F s ({q a J},{€"}) 



H(abaj3, cdyS) 



(29) 

90,91 



d €b d( fcd 

The first argument of Ti (4 components (aba/3) when a ^ b and 3 components (aa/3) 
when a = b but a ^ 0) denotes the index of the row of the matrix; the last one the 
column index. Because of the symmetry of the order parameters (^0) we can always 
take a < b or a < (3 when a = b. Therefore the square matrix 7i has dimension 
|[rn(n — 1) + r(r — l)ri 2 ]. One can distinguish three groups of matrix elements: firstly, 
those related to RSB fluctuations around q\ only, 



A, = H(aa(3,aa(3) = - J + J 2 ■ Sf) 2 ^ - q 2 

A 2 = H(aa(3, aa6) = H(aa(3, aj(3) = J 2 {((S£ • Sf S a a ■ sfy - q 2 
A 3 = H(aa(3 } ayd) = J 2 {((S£ • Sf sf 7 ^)) ~ q 2 ) 

A 4 = H(aa(3 } C1 5) = J 2 • S* SI ■ S s c )} - g 2 } ; (30) 

secondly, those related to fluctuations around go, 



B 1 = H{aba(3, aba(3) = -B + B 2 \({(S a a ■ Sff\ ) - q\ 



B 2 = H(abaP, abaS) = H{aba(3, abyP) = B <((S° ■ S b S%- S b )) —qf 



B 3 = H(abaf3 } abyS) = B 2 \((S a a - S p b ST a ■ S s b ) ) - q\ 



B A = H(aba(3 } adad) = H(aba(3 } dry/3) = B <((S° ■ S b ■ S d )) — qf 
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B h = H{aba/3, ad^S) = H{aba(3, c6 7 5) = B 2 |7(S£ • S% SI ■ S s d )^ - gjj 

B e = H(aba(3, cd^5) = B 2 {((s a a ■ S% sf^sfj) - qf) ; (31) 
and finally the matrix elements describing mixed RSB fluctuations 



d = H(aa0, ada5) = H(aap, cajP) = JB M(S" ■ Sj S% ■ S S d ) J - q qi 

C 2 = H(aaf3, adj5) = H(aaf3, ca 7 5) = J B M(S% ■ S% ■ S s d )^ - g gi} 

C 3 = H(aa(3 } cd 7 6) = JB ^(s a a ■ 5f ST e ■ S 8 d )} - q qi] , (32) 
with 

B = \B(3 2 J = \(B + J)(3 2 . (33) 

A simple interpretation of these matrix elements, similar to that given in e.g. |§, is not 
possible here, due to the vector character of the spins. 

The RS solutions are stable when the matrix (|29|) is negative definite. Upon 
analysing all eigenvalues (see [Appendix A|) it turns out that only two of these can 



cause the occurrence of a region of broken stability. The first replicon eigenvalue, which 
we will call the Greek replicon, reads 

X G = A 1 - 2A 2 + A 3 (34) 

and determines the Greek de Almeida- Thouless (AT) line Ac = 0. This AT- line measures 
the breaking of the symmetry with respect to the Greek indices. The corresponding 
eigenvectors are given by eq. (|A.11| ). The structure of the Greek replicon resembles the 



one of the replicon mode of the SK model found in |I7j], and also the Greek replicon 
mode of the SK model with coupled dynamics as studied in H. Since the Greek replicas 
are (by construction) related to the spin dynamics, the associated region of broken 
symmetry is located in the region of the phase diagram with low temperature T. The 
Roman replicon eigenvalue is given by 

X R = (B 1 - 2B 2 + B 3 ) + 2n(B 2 -B 3 -B 4 + B s ) + n 2 (B 3 - 2B 5 + B 6 ) (35) 

and measures the breaking of the symmetry with respect to the Roman replicas. 
This occurs at low coupling temperature T. Similar to the Greek replicon (0), the 
eigenvectors corresponding to the Roman replicon instability are symmetric under 
interchanging all but exactly two - in this case Roman - indices. It turns out that 
eigenvalues corresponding to eigenvectors which are symmetric under interchanging all 
but a number of indices which is larger than two, whether Roman, Greek or a mixture 
of both, can not induce an extra region of broken replica symmetry. Therefore, all 
regions where the RS Ansatz is unstable, are defined by the locations of the Greek 
AT line (Ac = 0) and the Roman AT line (Xr = 0). These lines are drawn in Fig. [I] 
and Fig. |3|. The latter shows explicitly that there is no re-entrance from the region SGI 
to the region with broken replica symmetry (RSB) when T is varied for fixed T. The 
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Figure 3. Phase diagram of the XY spin glass with slow dynamic couplings, drawn 
in the T-T plane; for Bq = 0, B = 1 and J = 3. Further notation as in Fig. [j]. 



RS solution is always stable in SGI with respect to the Roman replicas. In fact one can 
show analytically that the Roman AT line coincides with the SG1-SG2 transition line. 

The RS replica theory developed for our model, with spin and coupling dynamics 
on two different time-scales, is reminiscent of that of the simple XY model with one step 
replica symmetry breaking (1RSB). Our eigenvalues also formally resemble e.g. those 
describing the stability of the 1RSB solution in the perceptron model JTI2], ^3]. Note 
also that the position of the Roman AT line in our model is quite different from that 
in P|, although the phase diagrams of both models are qualitatively the same. The set 
of eigenvectors given there turn out to satisfy only part of the required orthogonality 
relations used in their calculations. An improved phase diagram for the SK model can 
be found in [|l^ . 

Finally we remark that the simpler model with homogeneous biases, mentioned 
earlier, does not involve Roman replicas, such that there appears only a Greek AT line. 
The latter line is qualitatively the same as the one in the model considered here. 



5. Simulations 



In order to complete our study and verify the predictions of our theory, we have 
performed numerical simulations of our model (note that, due to the parameter 
redundancy in our model, we can always restrict ourselves to J = 3 and B — 1). We 
have considered a population of XY spins, evolving according to the coupled Langevin 
equations (^|) and (|), which were discretised according to a standard Euler method, 
with iteration time step At = 0.001. Several interesting and subtle aspects arise when 
one attempts to carry out numerical simulations of models of the type studied here, 
with its widely disparate time-scales. Firstly, it will be clear that the presence of 
two adiabatically separated time-scales induce extremely large computing times, which 
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Figure 4. Evolution in time of the configurational energy (g) of the system, for 
parameters Bq = 0, B = 1, J = 3, T = 1.1 and n = 5, and with a system of size 
N = 200. The first window is chosen on the basis of the time required for the combined 
dynamical system (spins and interactions) to reach equilibrium; here we decided on 
a window size of 200. The values of the observables qo and q\ were obtained by 
performing a temporal average over a (third) time-window of the same size. 

prevent us from numerical exploration of the equilibrium regime for large system sizes. 
This is a general and systematic constraint, which causes important finite size effects, 
mainly near the phase transitions. In all our numerical studies we have, as a result, 
been forced to restrict ourselves to relatively modest systems of N = 200 spins. A 
second point concerns the evolution of the relevant quantities of the problem, spins and 
couplings, in view of the need to calculate the two main observables of the problem 
through an averaging process. This will have to be done very carefully, in order for the 
measured objects to indeed be identical to (or at least an acceptable approximation of) 
those calculated in the theory. Again the problem is related to having a finite system 
size: this narrows the time window where, on the one hand, the fast processes can 
be assumed to have been equilibrated, yet, on the other hand, the slow processes can 
be assumed not to have taken place. Thirdly, there is the fundamental problem that 
in regions where replica symmetry no longer holds (beyond either of the two AT lines) 
already the spin dynamics will exhibit the traditional phenomena associated with ageing, 
including extremely slow relaxation towards equilibrium; even without the additional 
superimposed slow dynamics of the couplings, it would have been extremely difficult to 
carry out numerical simulations that would probe the true equilibrium regime. 

We have dealt with these practical problems by adopting the following strategy. 
For a given set of couplings we first let spins to relax to their stationary state; then we 
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Figure 5. Spin-glass order parameters qo (circles) and q\ (squares) versus 
temperature, for n = 2. Continuous lines represent the theoretical predictions, while 
symbols denote simulation results (with N = 200, averaged over the time- window 
indicated in figure || and over 10 samples). As in the previous figures: Bq — 0, B = 1 
and J = 3. 
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Figure 6. Spin-glass order parameters qo (circles) and q\ (squares) versus 
temperature, for n — 5. Continuous lines represent the theoretical predictions, while 
symbols denote simulation results (with N = 200, averaged over the time- window 
indicated in figure || and over 10 samples). As in the previous figures: Bq = 0, B = 1 
and J = 3. 
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Figure 7. Non-normalised phase distribution P{<j>) = J2i ^[0 ~ 4>i] (discretised to 
a histogram) as observed at three different stages of the dynamical process towards 
equilibrium, for parameters Bq = 0, B = 1, J = 3, T = 1.1 and n = 5, and with a 
system of size N = 200. Bottom graph: (random) phase distribution at t = 0. Middle 
graph: phase distribution at t = 100 (during transient stage). Top graph: phase 
distribution at t = 200. 

perform the average (Si ■ Sj) over a number of time steps sufficiently large to have a 
statistically reliable measurement. We subsequently modify the interactions {Jij} for a 
certain number of time steps, completing what we call a 'dual iteration step'. This dual 
process is repeated until the global equilibrium state is reached. The key questions in 
the adequate employment of this strategy is to quantify rationally the various durations. 
According to the theory, since the time-scale associated with the couplings is infinitely 
slow compared to that of the fast variables (the spins), in each dual updating step we 
should modify the interactions {Jij} only very slightly (and during only a small number 
of update steps At). However, there are limits in practice to the extent to which 
one can proceed in this manner, in view of the danger of the simulations becoming 
so slow that they exceed by far one's computing resources. In our simulations we have 
updated the interactions during several hundred steps At (after having satisfied ourselves 
experimentally that a duration somewhere between 500 and 1000 iteration steps is quite 
appropriate) before, in turn, allowing the spin states to evolve. In this manner we 
have managed to speed up the convergence process towards global equilibrium, whilst 
continually verifying that the stationary values of the order parameters go an d qi, thus 
obtained, are not significantly affected. Even more delicate is deciding on the amount 
of time during which to evaluate the spin averages occurring in the stochastic equations 
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for the interactions. If the number of time steps used to calculate these averages is 
too large, the spins will have enough time to diffuse over the whole circle (due to finite 
size fluctuations which would have been absent in an infinitely large system), leading 
to an underestimation of go- Our experiments indicate that averaging over a period 
of between 2000 and 3000 iteration steps (of duration At each) gives reliable results. 
Finally, we have to decide on the window size (the number of dual updating steps) 
which we have to average over in order to compute the observables of the system. The 
logical approach would appear to be to monitor the evolution of quantities such as the 
energy (fj), starting from the initial state, until the stationary state has been (or at least 
appears to have been) reached. Figure |4| shows a typical numerical experiment. The 
dynamics towards equilibrium on this time-scale is ultimately controlled by the slow 
variables, the couplings. The spins respond to changes in the couplings in a stochastic 
master/slave fashion, and only when the slow variables (the couplings) have reached a 
stationary stationary state can we speak about global (thermal) equilibrium. In figure 
£| we see that 200 dual steps suffice to ensure the absence of the main transient effects. 

The observed dependence on temperature of the two spin-glass order parameters, 
go and q±, is illustrated in figures § and || together with the corresponding theoretical 
predictions. We have carried out these numerical simulations for the temperature 
ratios n = T/T = 2 and 5, respectively. For the smaller temperature ratio n = 2 
we observe that our simulations indeed confirm the existence of two spin-glass phases; 
one exhibiting freezing only on spin time-scales, and a second spin-glass phase where 
freezing is observed on all time scales. However, quantitative agreement between theory 
and simulations is extremely difficult to achieve, due to the practical problems outlined 
above. For n = 5 the system is in the region where the theory predicts that a first 
order phase transition from a paramagnetic phase to a second-spin glass phase should 
be found. We observe a good agreement between theory and simulations, except for 
temperatures close to the transition, where finite size effects are obviously increasingly 
important. In addition to the above equilibrium observables, we have investigated 
other, non-equilibrium, aspects of our model, by way of further illustration. We have 
measured, for instance, the distribution P{4>) = J2i S[4> ~ 4>i\ °f phases (pi, defined via 
Si = (cos 9i, sin 9i). Figure [7] shows this distribution at three different stages during 
the evolution towards equilibrium, for system parameters identical to those of figure 
Initially the phases were distributed uniformly; one observes this distribution to deform 
spontaneously into a bi-modal one, driven in conjunction with the feedback provided by 
the (slow) dynamics of the couplings. 

6. Concluding Discussion 

In this paper we have discussed and solved a version of the classical XY spin-glass 
model in which both the spins and their couplings evolve stochastically, according to 
coupled equations, but on widely disparate time-scales. The spins play the role of 
fast variables, whereas the couplings evolve only very slowly, but according to local 
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stochastic laws which involve the states of the spins. In the context of disordered 
magnetic systems this model describes a situation where one takes into account the 
possible effects of slow diffusion of the magnetic impurities, without necessarily assuming 
energy equi-partitioning between the slow variables (the impurity locations) and the 
spins (hence the potentially different temperatures associated with each). Alternatively, 
in the context of neural systems this model would describe coupled neural oscillators 



15| with autonomous stochastic Hebbian-type synaptic adaptation on the basis of the 



degree of firing synchrony of pairs of neurons. 

We have solved our model within the replica-symmetric (RS) mean-field theory, 
involving two levels of replicas: one level related to the (slow) couplings, and one 
related to the disorder in the problem (the symmetry breaking terms in the dynamics 
of the couplings, representing preferred random values of the latter). The solution of 
our model, in RS ansatz, is mathematically similar to that of the XY model with static 
couplings, but with one-step RSB. This is reminiscent of the general connection between 



the breaking of replica symmetry and the existence of dynamics on many time scales [24 



We have discussed in detail the stability of the RS-solutions, including the calculation of 
all eigenvalues and their multiplicities (details of which can be found in the appendices). 
It turns out that two distinct replicon eigenvalues determine the region of stability, and 
thus the region of validity of the RS solution. 

The thermodynamic phase diagram is found to exhibit two different spin-glass 
phases, one where freezing occurs on all time-scales, and one where freezing occurs 
only on the (fast) time-scale of the spin dynamics. We also find both first- and second- 
order transitions; the origin of the first order ones is the positive feedback in the system 
(compared to a system with stationary spin-couplings) which is induced by the super- 
imposed coupling dynamics. As could have been expected, the physics of the present 
model resembles that of the SK model with dynamic couplings, apart from a re-scaling in 
temperature and provided an appropriate adjustment of the calculation of the AT lines 
in || is made. Our calculations show how the methods used for solving the Ising case can 
be easily adapted to deal with more complicated spin types, and in addition illustrates 
further the robustness of the phase diagrams describing the behaviour of large spin 
systems with dynamic couplings. 

Numerical simulations present further interesting technical challenges, due to the 
existence of adiabatically separated time scales (which requires equilibration of two 
different nested stochastic processes in order to test the theory), in addition to the 
already highly non-trivial and extremely slow dynamics of the fast (spin) system. In spite 
of the important finite size effects, which are inevitable given the practical constraints 
on available CPU time, our results show good agreement with the theory and confirm 
the main characteristics of the predicted behaviour. 
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Figure Al. Graphical representation of the structure of a general eigenvector of 
the Hessian matrix (|2^). Dark spaces denote the e-components; the other non-empty 
spaces are the 77-components. Further details are found in the text. 



Appendix A. The Hessian matrix 

Appendix A. 1. Eigenvectors and eigenvalues 

In this Appendix we show how to find all eigenvectors, eigenvalues and their multiplicity 
of the Hessian matrix (p9|). We immediately remark that 7i is symmetric, implying 
that its eigenvectors must be orthogonal, a property exploited heavily in finding its 
eigenvalues. 

We denote the eigenvectors in the ^[rn(n — 1) + r(r — l)n 2 ]-dimensional space by 
(e"' 3 , rfcd) with a < (3 and c < d, and represent them graphically in a square matrix 
as in Fig. [Al| . This matrix is of size nr x nr and is divided in r 2 sub-matrices of size 
n x n, labeled by two Roman indices (a,b,c = l...r). The elements of these sub- 
matrices, in turn, are labeled by two Greek indices (a, f3,j,5— 1 . . . n). Thus the rows 
and columns of the matrix carry a Roman and a Greek index. Each of the matrix 
elements corresponds to a component of the vector {e^f, 77^), except for the diagonal 
components, which are put equal to (viz. an empty space). The elements of the type 
(act, cry) correspond to r]^ when a 7^ c and to e" 7 when a = c but a 7^ 7. The matrix is 
symmetric, such that rf^ = v2a an d e" 7 = e cT- 

The eigenvectors of 7i with eigenvalue A satisfy the eigenvalue equation 

w (S) =A (l)' (A - i) 

At this point we remark that {e^ 13 } are generally uncorrelated fluctuations of the order 
parameters {q^ 13 } and {rfl^} of {q2d}- Therefore we call a vector (e^r/^) symmetric 
under permutation of the Roman and Greek indices when the components {e"^} and 



Coupled dynamics in the XY spin-glass 



20 




Figure A2. Representation of the replica symmetric eigenvector. Empty space 
denotes zero elements; spaces with the same fill pattern denote elements with identical 
matrix elements. The dark spaces indicate the e-components; the other spaces indicate 
the ^-components. 



{rfl^} are simultaneously symmetric under permutation of these indices. In order to 
find the explicit form of the eigenvectors we make a general proposal based on this 
symmetry. Furthermore, we use the eigenvalue equation (|A.1|) and the orthogonality of 
the eigenvectors corresponding to different eigenvalues. 

We start from the symmetric solution where all components are identical, viz. 

ef = f and rj£ = g . (A.2) 

These vectors are represented in Fig. [A"2| . Substitution of ( |A.2| ) into the eigenvalue 
equation ( A.l|) reduces the number of equations to solve to two and we easily find 

(X 1 - A) / + Yi g = 
X 2 f + (Y 2 - A) g = 

Xi = Ai + 2(n - 2)A 2 + -(n - 2)(n - 3)A 3 

Y 1 = 2n(r - l)d + n{n - 2)(r - l)C 2 + Ki 2 (r - l)(r - 2)C 3 . (A.3) 

From this we get two non-degenerate eigenvalues 
1 



Ai, 2 = -[Y 2 + X 1 ± y/(Y 2 + X,y - 4(X 1 y 2 - X 2 Y 1 )j (A.4) 

X 2 = 2{n - l)C x + ( n - l)(n - 2)C 2 + Ki{n - l)(r - 2)C 3 
Y 2 = B x + 2(n - l)S a + (n - l) 2 fi 3 + 2n(r - 2)B 4 

+ 2n(n - l)(r - 2)5 5 + ^n 2 (r - 2)(r - 3)^6 . 
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Figure A3. Representation of the eigenvectors ( A. 6 ). Empty space denotes zero 
elements; spaces with the same fill pattern denote elements with identical matrix 
elements. The dark spaces denote the components e^P , the striped spaces the e"' 3 , 
the checked spaces the rfj^, the other non-empty spaces the T]"^ . 



The matrix elements A±, . . . , C3 are given by eqs. 

The other eigenvalues are related to the breaking of the symmetry in ( |A.2|) , both 
at the level of Roman indices and at that of the Greek indices. The most simple form 
of symmetry breaking is the case where almost all components are identical, except for 
those labeled by a single specific pair of indices {x, 9}: 

Vxb = k ; Vxb = l\ Vab = m a , b + x; a, (5 ^ 9 (A.5) 



This increases drastically the number of equations obtained from ( |A.1| ). A first trial 



solution within this group of candidate eigenvectors is obtained upon putting h = m = 0, 
giving a group of eigenvectors with all components vanishing, except the ones labeled 
by {x,9}: 

6(3 _ f . a/3 = _ 2 f . _ n . 

c x J 1 c x ^ 2 ^ ' a ' 

Vxb = k i Vxb = 7 ^ ) Vab = ® 

n — 1 

Y x k = {X-Xl) f 

X 1 = A 1 + (n- A)A 2 - (n - 3)A 3 

Yi = (n - 2)(r - 1) (d - C 2 ) . (A.6) 

n — 1 

The graphical representation of these eigenvectors in the form of a matrix is drawn in 
Fig. [A3[ The associated eigenvalues read 



A 3 ,4 = - ( Y 2 + Xt ± y/(Y 2 + Xi) 2 - A{X X Y 2 - X 2 Y\ 
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Figure A4. Representation of the eigenvectors ( A. 8 ). Empty space denotes zero 
elements; spaces with the same fill pattern denote elements with identical matrix 
elements. Striped spaces denote the components t^P , dark spaces the components 
, checked spaces the components 7]"^ , the other non-empty spaces the components 

ad 
Vab ■ 



X 2 =(n-l)(d-C 2 ) 

Y 2 =B 1 + (n- 2)B 2 - (n - 1)B 3 + n(r - 2)B 4 

+ n(n-l)(r-2)B 5 . (A.7) 

The degeneracy of the associated eigenspace is r(n — 1); it is found by calculating 
explicitly the rank of the matrix composed by these eigenvectors, as will be outlined in 
appendix A2. 

Insertion of the general proposal ( |A.5| ) into eq. ( |A.1| ), and subsequently requiring 
the orthogonality of this vector to the eigenvectors we have already found earlier, leads 
us to the new eigenspace 

1 



M _ 


f af3 _ 
c x 


/; 


c a 


Vxb = 


a/3 
Vxb — 


fc; 


a/3 
Vab 


Y l k -- 


= (A- 


XI) 


f 



r — 1 
2 

~r -2 



/; 

k 



X X = A X + 2(n - 2)A 2 + - (n - 2)(n - 3)A 3 - -n{n - 1)A 4 

Yi = 2n(r - l)d + n(n - 2)(r - 2)C 2 - n 2 (r - 1)C 3 (A.8) 



with eigenvalue 



1 



A 5 , 6 = - ^ 2 + x x ± ,j(y 2 + x x f - 4(x 1 y 2 - x 2 Yi 

X 2 = ^| f (n - l)d + ^ (n - l)(n - 2)C 2 -\n{n- 1) C 3 ' 
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Figure A5. Representation of the eigenvectors ( A. 11 ). Empty space denotes zero 
elements; spaces with the same fill pattern denote elements with identical matrix 
elements, the striped spaces indicate the components e£", the dark spaces the 
components e e x a and e" u , the checked spaces the components t^P . 



Y 2 = B X + 2(n - l)B 2 + (n - 1) 2 5 3 + n(r - A)B. 



i 



+ n(n-l)(r-A)B 5 -n 2 (r-3)B 6 . (A.9) 

The vectors ( |A.8| ), represented graphically in Fig. |A4| , are symmetric under interchanging 
all indices but one Roman index x, and have associated with each eigenvalue a (r — 1)- 
dimensional eigenspace. 

A second class of eigenvectors is found by considering a situation where one Roman 
index x and two different indices 6 and v cause breaking of the replica symmetry. In 
their most general form these are given by 

e ic / ) ^-x ^x 9 i ^x^ ^ ) k ) 

Vxb=Vxb= l 'i Vxb=m; Vab=P {a,b ^ x; a, (3 ^ 6,v) . (A.10) 

More explicitly, we propose a vector with two special Greek indices, i.e. we try a solution 
with k = l = m = p = 0. It corresponds to a vector with all 77-components and all 
e-components which are not related to this Roman index vanishing, and with broken 
replica symmetry with respect to the two Greek indices: 

1 2 

0u _ r . 0(3 _ a,v f . a/3 _ r . a/3 _ n. 

x 1 ' x x " n-2 1, x ( n -2)(n-3) 1, a ' 
ril = r$ = 0\ Cf = 0; rtf = (a,b^x; a,p^6,v) . (A.ll) 
They are visualised in Fig. ~KE . The eigenvalue equals 

X 7 = A 1 -2A 2 + A 3 (A. 12) 

and has multiplicity rn(n — 3)/2. 
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Figure A6. Representation of the eigenvectors ( A. 15 ). Empty space denotes zero 
elements; spaces with the same fill pattern denote elements with identical matrix 
elements. Dark spaces denote the components 77"^ , striped spaces the components 
r]™P and 77"^, the other non-empty spaces the components 77"^ . 



Finally, replica symmetry can also be broken by two spins with different Roman 
indices x and y. Upon calling the corresponding Greek indices 9 and u, this group of 
eigenvectors reads in its most general form 



rW3 



aB 



c aB 



Vxy ~ k i Vxy ~ Vxy 

aB olB aB ol8 

Vxy = P ; Vxb = Vay = Q \ Vab 



68 



y -9; 



c aB 



h; 



m ; 



Vxb — Vay 

t (a,b 7^ x, y; a, (3 7^ 9, u) . (A. 13) 

Again ( |A.1|) and the orthogonality relations are used in order to find explicit solutions. 
One eigenvalue is given by 

A 8 = (B 1 - 2B 2 + B 3 ) + 2n(B 2 -B 3 -B 4 + B 5 ) + n 2 (B 3 - 2B 5 + B 6 ) (A.14) 
with the corresponding eigenvectors (Fig. |A6|) 



JB = u8 = aB 
c x c y c 



88 _ av _ aB _ aB 
Ixb lay Ixb lay 



6v 



- a P = f a @ = D ■ ri"~ = n' r ' = n~~ = n 

r '-y c a w ) Ixy Ixy Ixy Ixy 

1 



aB 



k: 



aB 
Vab 



r-2"' ,ao (r-2)(r-3) 

and with degeneracy r(r — 3)/2. All e-components vanish. 
Yet another eigenvalue is 

A9 = B\ — 2B 2 + B 3 

with eigenvectors given by 



(A.15) 



(A.16) 



V%) — k ] 7], 



V, 



n 



V. 



aB 
xy 



(n-lf 



k: 



Coupled dynamics in the XY spin-glass 



25 



v 





Figure A7. Representation of the eigenvectors ( A. 17 ). Empty space denotes zero 
elements; spaces with the same fill pattern denote elements with identical matrix 
elements. The dark spaces denote the components if" the striped spaces the 



components rfj^ and r]"y, the other non-empty spaces the components rj" 



a/3 



c"0 



= 

c x 

aB 
Vxb 



c aB 



Vxb = Vay = Vxb = Vay = Vab = 5 ( A - 17 ) 

This corresponds to a situation where only the components with the marked indices x 



c a/3 _ 
-a 

aB 
: Vab 







and y are non-vanishing. The vectors are drawn in Fig. |A7| . The degeneracy of the 
eigenvalue (|A.16| ) is \ r(r — l)(n — l) 2 . 



Finally, the last eigenvector reads (Fig. [AS 



n 9v 

Ixy 



7i aP 
'Ixy 



y 

k ■ n 913 

™ ) 'Ixy 
1 



e af3 = e af3 



-y 

av 
Ixy 



n 







aB 
Vab 



. 



2(n 



'Ixb — 'lay 



2(r - 2) 



n 



r 



k', Vxb 



ay 2(n-l)(r-2) 



The corresponding eigenvalue is equal to 

A10 = B x + (n - 2)B 2 - (n - l)B 3 - nB A + nB 5 
and has degeneracy |r(r — 2)2(n — 1). 



(A.18) 



(A.19) 



Appendix A. 2. The multiplicity of the eigenvalues 

This appendix has been included in the present paper since no information is available 
in the replica literature about explicit methods to find the multiplicity of the eigenvalues 
of the Hessian matrix. We do not aim for mathematical rigour, but just aim to aid the 
reader by giving a heuristic method for finding the solution. 



Coupled dynamics in the XY spin-glass 



26 




Figure A8. Representation of the eigenvectors ( A. 18 ). Empty space denotes zero 
elements; spaces with the same fill pattern denote elements with identical matrix 
elements. The dark spaces denote the components ff^L, the checked spaces the 
components if® and 77"^, the diagonally striped spaces the components 77"^, the 



8 S 

vertically striped spaces the components rjj^ or 77"^, the other non-empty spaces the 



components 77"^ or 77"^ 



-a/3 



Appendix A. 2.1. Eigenvalues X^^ and Xf,$ We first focus on As^. The sum of all r 
eigenvectors (|A.8| ) equals zero, such that there are at most r — 1 linearly independent 
vectors. In the sequel we show that the rank of the matrix composed of all eigenvectors is 
exactly given by this latter value. We consider the r-dimensional sub-matrix constructed 
by that part of the e-component of the vectors (|A.8|) with fixed Greek indices (e.g. 
a = 1, (3 = 2) Viz. (e^ 2 , a = 1, . . . , r); the r vectors are obtained by varying the Roman 
index x. This matrix reads 

M = (f-g) 1 + gJP g = ^-f (A.20) 

r — 1 

where 11 and IP are the unit matrix and the projector matrix, respectively, here both of 
dimension r. Since these matrices commute they can be diagonalized simultaneously. 
The only eigenvector of IP with a non-zero eigenvalue is the vector (1, 1, ... , 1) with 
eigenvalue r. All other eigenvectors can be chosen orthogonal to this vector. This results 
in the following eigenvalues for M: one non-degenerated eigenvalue (/ — g) + rg = 
and an (r — l)-fold degenerate eigenvalue / — g = ~zrf- We conclude that all four 
eigenvalues (|A.9|) have degeneracy (r — 1). 

The multiplicity of X^^ can be found in an analogous way. Since eigenvectors with 
a different Roman index x are always independent, it is sufficient to determine the 
dimension of the matrix constructed from the n eigenvectors with the Roman index 
fixed, e.g. x = 1, and with the Greek index 9 running from 1 to n. 
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Appendix A. 2. 2. Eigenvalues A 7 and Ag We start with the less complicated calculation 
for Ag. As before we consider a sub-matrix which is the r(r — l)/2-dimensional matrix 
of the 7/-part of the eigenvectors ( |A.15 ), with fixed Greek indices e.g. a = 1 = (3 viz. 
, a<b=l,...,r) 

M = (f-g)t + gW + (h-g)lB (A.21) 
-1 f 2 
9 r-2 1 (r-2)(r-3) 

B Q/3 , 7 5 = (1 - <L 7 )(1 - S aS )(l - 6fa)(l - 5p S ), a < (3,7 <5. 

Because the three matrices appearing above commute, the problem of finding the 
eigenvalues fii of M reduces to finding the eigenvalues 7, of IB. First, one finds the non- 
trivial eigenvector of IP, (1, 1, . . . , 1), which is also an eigenvector of IB with eigenvalue 
(r — l)(r — 3)/2. The other eigenvectors (x a b), with a < b, can be chosen orthogonal to 
the trivial one, i.e. J2 a <b x ab = 0, leading to the following simplified eigenvalue equation 
for the eigenvectors of IB 

1 x ab = ^ab,cd x c d = x a + y a + x b + y b + x ab 

c<d 

JJ- X a = Y X ab y b = Xab 

b(>a) a(<b) 

(1 - 7)x afe = x a + y a + x b + y b . (A.22) 

The first solution of these equations is 7 = 1, with multiplicity r(r— l)/2 — r due to the 
condition x a + y a = 0, a = l,...,r. When 7 7^ 1, we can sum (|A.22|) in two different 

ways: 

■ (1 - l)Vb = (6 - l)(ar* + 3ft) + XI ( x « + 

a(<6) a(<6) 

X : (1 - 7)x a = (n - a)(x a +y a ) + ^ (x 6 + y b ) ■ 

b(>a) b(>a) 

Adding the two equations gives 

(3- 7 -r)(x a + y a ) = (A.23) 

and leads to the eigenvalue 7 = 3 — r. Because of the second factor in (|A.23[) we can 
conclude that IB has no other eigenvalues than the ones already found, leading to the 
multiplicity r — 1 for 7 = 3— r. It turns out that only 7 = 1 gives a non-zero eigenvalue of 
M. We can therefore conclude that the rank of M is equal to r(r — 1)/2 — r = r(r — 3)/2. 

The same procedure can be followed to determine the multiplicity of A7, but now 
by considering M as composed by the e-part of the eigenvectors, for a fixed choice 
of the Roman index x. Eigenvectors with different Roman indices are again linearly 
independent. 



Appendix A. 2. 3. Eigenvalue A9 First, we note that each choice of the Roman indices 
x and y gives a set of independent eigenvectors, so we can limit ourselves to finding the 
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rank of the matrix generated by the vectors with x = 1 = y 
, ! " ' v/i=Vn=9 riS = h 



n-\ 

which can be written as 



< = / 



(n - If 



0)B 



a = 7 

a ^ 7 with B = TP 



(A.24) 



M = (/-</) l + gW + (h 
Bq/ ^ = { %s 

where the dimension of B is n. The special structure of B allows us to find its eigenvalues 
quite easily, namely: one non-degenerate eigenvalue (n— l) 2 , one 2{n— l)-fold degenerate 
eigenvalue — (n — 1), and one (n — l) 2 -fold degenerate eigenvalue 1. The first and second 
of these give a zero eigenvalue for M, leading to rank M = (n — l) 2 . 



Appendix A. 2. 4- Eigenvalue X 10 Calculating the multiplicity of this eigenvalue involves 
some more work. In contrast to the previous sections, we here have to calculate the 
dimension of the full matrix of eigenvectors, rather than just the dimension of a suitable 
sub-matrix. Upon writing the eigenvectors as 



'X 

rfi"-. 

Ixy 
'Ixy 



9 



y 

f; 

k; 

nr-2 



r aj3 



77 

'Ixy 

a/3 
Vxb 



f 



y 

'Ixy 
'lay 

h- 



9 



= ; 

6/3 
Vxb - 



a/3 
Vab 



lay 



. 

h; 



1 



/ 



2(n-l) " 2(r-2 / 
the matrix of eigenvectors M reads 

/ 
9 
k 
h 
I 



k 



(n-1) 



f 



2(n-l)(r-2) 



/ 



M 



abaf3,cd-yS 







a = c and b = d and a = 7 and (3 — 8 
a = c and b = d and (a = 7 or (3 = 5) 
a = c and b — d and a 7^ 7 and (3^5 
one subindex and the corresponding superindex are equal 
one subindex is equal (and not the superindex) 
otherwise 



(A.25) 



It is found to consist of sub-matrices B (containing elements /, g and k) on the diagonal, 
and sub-matrices HDj, i — 1, ... ,4 (containing elements h and /, ordered in four different 
ways) elsewhere. All of these sub-matrices have dimension n 2 . 

In order to simplify the problem we first construct the matrix C which is made up of 
columns which are orthogonal and normalized eigenvectors of B. For every sub-matrix 
IDj we construct C IDj C T , where C T is the transposed matrix of C. Next we construct 
the matrix 

/ C \ 



C 



(A.26) 



C ) 
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Since CC T = 1 by construction, the eigenvalue equations of M and M. = C JM C T are 
the same, and we can restrict ourselves to solving the simpler eigenvalue problem of the 
latter matrix. 

First we focus on the matrix IB and the construction of the matrix C . As in 
Section [Appendix A.2.3| we have 

JB = (f-g) l + gW + (k-g)JK 

JK = I ° a = 7 

a ^ S \ WLps a ^ 7 with K = P - 1 • 

Upon using the Gramm-Schmidt procedure for constructing a set of orthogonal and 
normalized eigenvectors of K, we arrive at the following result: 



C 



i 

1 

V2 

1 fl 

2 V 3 

I fl 

3 V 4 



1 
\/n 
1_ 

V2 

- 

2 V 3 

I /I 

3 V 4 



1 





\ n-1 



n-1 



1 

n— 1 



n-1 



1 

n— 1 



1 






1 







n-1 



1 

n— 1 



n-1 



1 

n— 1 



n-1 



1 







\ 



n-1 



(A.27) 



Due to the similar structure of the two matrices IK and IK, we now can immediately 
read off the matrix of eigenvectors C of IB: one takes a matrix with the structure of 
( A.27 ), and multiplies each matrix element by C, arriving at a matrix with dimension 
n 2 . Given this matrix it is straightforward to calculate C IDj C T for all sub-matrices 



ID,. We arrive at 



M 



abaft, cdjS 



2(n-l) 



2(n-l)(r-2) 



a = c and b = d and a 
a = c and b = d and a 
a = c and b ^ d and a 
a = d and b ^ c and a 







a/d and b 
a c and b ■ 
otherwise 



c and a 
d and a 



7 = 1 and j3 
7 7^ 1 and (3 
7^1 and (3 
6^1 and (3 ■ 
5=1 and (3 ■ 
7 = 1 and /5 



-5^1 

■ 5 = 1 
5 = 1 
7 = 1 

■ 5^1 



(A.28) 



In view of the large number of zero-rows in this matrix, it is convenient to define a 
|r(r — l)2(n — l)-dimensional matrix .M, which contains only the non-trivial rows 

■^^■aba^cdp A^ ab"f&,cd^iu 

1 5 = a + 1 for a = 1, . . . , n — 1 

a — (n — 2) 5 = 1 for a = n, . 

1 i/ = /3 + lfor/3 = l,. 

(3 - (n - 2) i/ = 1 for /3 = n, . 



7 
7 



(A.29) 



This matrix can be written as A4 



,2(n-l) 

,71—1 

,2(n-l) 

IN, where IN is a matrix with 



2(n-l) 1 2(n-l) (r-2) 

elements or 1 only. As in Section Appendix A. 2. 2 the eigenvalues of this matrix are 
obtained by summing the eigenvalue equation for IN in two different ways. We then find 
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the eigenvalue A = — 1 with multiplicity |r(r — 2)2 (n — 1), and the eigenvalue A = r — 2 
with multiplicity r(n — 1). It turns out that only the first of these eigenvalues gives a 
non-zero eigenvalue for Xi, and we may conclude that rank M = r(r — 2)(n — 1). 



